Is the Concentration of Dark Matter Halos at Virialization 

Universal ? 

Massimo Ricotti 1 , Andrew Pontzen 2 , and Matteo Viel 2,3 
ABSTRACT 

Several recent studies suggest a correlation between dark matter halo mass 
and the shape of the density profile. We re-analyze simulations from Ricotti 
(2003) in which such a correlation was proposed. We use a standard analysis of 
the halo density profiles and compare the old simulations to new ones performed 
with Gadget2, including higher resolution runs. We confirm Ricotti's result that, 
at virialization, the central log slopes a, at 5%-10% of the virial radius are cor- 
related with the halo mass and that the halo concentration, c, is a universal 
constant. Our results do not contradict the majority of published papers: when 
using a split power law to fit the density profiles, due to the a — c degeneracy, 
the fits are consistent with halos having a universal shape with a = 1 or 1.5 and 
concentrations that depend on the mass, in agreement with results published 
elsewhere. 

Recently, several groups have found no evidence for convergence of the inner 
halo profile to a constant power law. The choice of a split power law parame- 
terization used in this letter is motivated by the need to compare our results to 
previous ones and is formally valid because we are not able to resolve regions 
where the slope of the fitting function reaches its asymptotic constant value. 
Using a non-parameterized technique, we also show that the density profiles of 
dwarf galaxies at z ~ 10 have a log slope shallower than 0.5 within 5% of the 
virial radius. 

Subject headings: Galaxies: formation; Methods: N-body simulations; Cosmol- 
ogy: theory. 
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Introduction 



It is widely believed that dark matter halos have a universal density profile. In the last 
decade the debate focused largely on whether the slope of the density profile, a, approaches 
values close to 1 or 1.5 at small radii. The concentration, c, of dark halos supposedly 
incorporates all the physics: the dependence on the cosmological model, the halo mass and 
redshift of virialization. However, the determination of the inner slope of halos is uncertain 
due to the scatter of halo shapes and a well known degeneracy between a and c when fitting 
the density profiles. Indeed, only a few published simulations have sufficiently high resolution 
to partially resolve this degeneracy and these simulations only cover a limited range of halo 
masses (typically Milky- Way type or cluster type halos) at redshift z = 0. 



In this paper, we re-analyze data from N-body simulations performed in iRicottil (120031 ). 
henceforth R03 , and compare them to a new set of simulations performed using Gadget2 
( jSpringell 120051 ) at the same and higher resolution than in R03. The R03 results are of par- 
ticular interest because, in R03, it was first found that the log slope of the inner parts (at 
10% of the virialization radius for just virialized halos) of the dark matter density profile, a, 
varies with the mass of the halo. This result may ease some tension between theory and ob- 
serva tions. Indeed, for low mass halos (< 10 9 M(T, ), which should correspond to dSph galaxies 



Ricotti fc GnedirJ 120051: iRead et al.l bopdT and perhaps some Low Surface Brightness 



[e.g. 

(LSB) galaxies (e.g., Ide Blok fc Bosmall2002l ). R03 finds a < 1 at 10% of the virial radius 
(where r V i r ~ 1 kpc). In R03, Milk y- Way type halos w ith mass M ~ 10 12 M & are well 
fitted by a NFW profile with a ~ 1 (INavarro et al.lll996l ) and there is also some evidence 
that high er mass systems (> 1O 15 M ) have steeper cusps, a > 1. This is consistent with 
results in iMoore et al.l (119991 ). although this latter paper also proposes a ~ 1.5 for rather 



small, galaxy mass halos (~ 10 M ). It is important to clarify that the quoted values of 
a in the inner part of the halos depend on the resolution of the simulations (e.g., in R03 is 
fixed to be about 10% of the virial radius). There is no reason to believe that a converges 
to any asymptotic value as suggest ed by the arguments in R03 and and by high-resolution 
simulation of Milky- Way size halos (INavarro et al.ll2004l ; iGraham et al.ll2005l ). Note that the 
result in R03 does not depend on this assumption. In R03 we show that NFW profile does 
not provide good fits for small mass halos at z ~ 10 by comparing the shapes of the circular 
velocities of just virialized halos with widely different masses (hence, virialized at different 
redshifts). Our conclusion is that density profiles do not have an universal shape as often 
assumed before. A theoretic a l inte rpretation of this result is also proposed based on previous 
work by ISubramanian et al.l (120001 ) in which a simple relationship that relates a to the slope 
of the power spectrum of initial density perturbations is provided. Our results are consis- 
tent with this simple scaling relationship, suggesting that the halo shape at a given mass or 
spatial scale depends on the the slope of the power spectrum at that scale. This result is 
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based on fo r mal fi tting; of the circular velocities of the halos with generalized NFW profiles 
Colin et al.l (120041 ) performed similar simulations and failed to find shallow core s in a similar 



mass range. Nevertheless, similar correlations have been found by oth ers, fe.g. iJing fc Suto 



2000 



and 



Taylor Sz Navarrol 



Merritt et al 



200ll ; ICen et all 120041 ) . More recent work bv iGraham et al.l (120051 ) 



(120051 ) also find a correlation between halo mass and the shape of the 
density profile. In this case they parameterize this correlation in terms of the Sersic index , 
n, rather than a central log-slope dependence. In (INavarro et al.ll2004l ; iMerritt et al.l 120051 ) 
is argued that this parameterization provide a better fit of high-resolution halos than a split 
power-law. 

The main motivation of this letter is to understand whether the R03 results are in 
contradiction with previous works and clarify if the discrepancy can be attributed to the 
method of analysis, the N-body code used by R03, or insufficient resolution of the simulation. 
Ricotti used a P 3 M N-body integrator and analyzed the data using circular velocities instead 
of density profiles. Here we adopt the widely used "tree" N-body integrator Gadget2, with 
256 3 and 512 3 particles and we develop a quantitative method to analyze the density profiles 
using a standard \ 2 minimization technique. 

This letter is organized as follows. In Section [21 we describe the set of simulations from 
which data have been used and the procedures we have adopted and developed for analyzing 
individual halos; in Section [3] we present the results; finally we conclude in Section HI 



2. Simulation Data and Analysis 



All simulations used and referred to in this work have identical cosmological parameters: 

0.91 and h = 0.7. R03 uses a P 3 M integrator 



Gnedin & Bertschineer 


1996 


i while 


( Springe! 


2005; 


Sprineel et al. 


2001 


)■ 



is: Run-L, where Run = R03, GR03, G256 and G512 describe a different simulation and 
L = 1,32,256 refers to the box size (i.e., 1, 32, 256/r x Mpc). The runs "R03" refer to the 
original R03 simulations, "GR03" use the initial conditions in R03 but are re-run using 
Gadget2, "G256" and "G512" are new runs using Gadget2 with 256 3 and 512 3 particles, 
respectively. The purpose of running the new simulations GR03, G256 and G512 was to 
ensure the results in R03 were not affected by any irregularity in the simulation method 
employed. GR03 checks for problems with the R03 simulation parameters and code; G256 
checks the initial conditions and G512 checks for resolution related issues. 



The redshift of analysis in all lh 1 Mpc simulations is z = 10. The larger box sizes are 
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analyzed when the clustering of the most massive halos is similar, which turns out to be z = 3 
for L = 32/i _1 Mpc and z = for L = 256/i _1 Mpc. We extracted the halos using a Friends-of- 
Friends algorithm with linking length I = 0.2 (chosen by analog y with the spherical collapse 
model). We follow the iterative method of iPorciani et al.l (120021 ) to ensure that the extracted 
halos are bound. At each stage the total energy of each particle is calculated. Particles which 
do not appear to be bound are excluded from the potential calculation in the next stage. 
An accurate determination of the halo centers is important, since from a miss-centered halo 
we would deduce a systematically flattened profile. We use three centering methods: i) 
Density Maximum. An algorithm for finding the density maximum uses an adaptive grid of 
cubes; where the number of particles inside a cube exceeds a certain threshold , the cube is 



subdiv ided and the process iterates, it) Shrinking Sphere Center of Mass, as in lPower et al. 



( 120031 ). Here we calculate the center of mass (COM) and then iterate, including at each stage 
only particles within some sphere around the previously calculated center. Hi) Potential 
Minimum. The particle with the minimum potential is chosen as the center. Generally we 
find that the potential minimum lies within the convergence radius of the shrinking sphere 
center. All three alternative centering procedures are used; if one procedure does not produce 
results within the convergence radius of each of the others, the halo is flagged as unusable 
for spherically averaged processing. We use the convergence criterion of iPower et al.l (120031 ) 
to determine a central region where the density distribution is unreliable. This yields a 
convergence radius inside which the profile is not considered during the fitting procedure. 

We fit th e density profiles usi ng a standard x 2 minimization routine with a generalized 
NFW profile (jNavarro et al.lll997l ) in which the concentration c and inner slope a are free 



parameters: p oc x a (x + 1) 



-3+a 



where x = r/r s and r s 



./c. It is obvious that 



this parameterization of the profile should be regarded valid only between the virial radius 
and the convergence radius which, for all the halos considered in this work, is a constant 
fraction ~ 5 — 10% of the virial radius. In this pa per, we use the Poisson error fo r our 
X 2 minimization procedure. It is widely noted (e.g. iJingi l2000l : iNavarro et al.l l2004f ) that 
errors introduced by deviations from an "idealized" equilibrium profile (i.e. substructure, 
asphericity and irregularities) are likely to be at least as important as Poisson errors. 



3. Results 

When fitting the split power law to a halo, there is a degeneracy between c and a: the 
model is effectively constrained not to a point but to a line through the c-a plane, along 
which x 2 varies by only a small factor ~ 2 from its minimum to the edges of the region of 
interest < a < 1.5. The location of locus of minima varies widely from halo to halo and 
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may select non-physical configurations with a < 0. Perpendicular to this line, however, the 
value of x 2 increases rapidly. In Figure [1] left panel, we illustrate this degeneracy for the 
10 most massive halos in R03-1. By constraining either c or a arbitrarily, each line can be 
reduced to a point, given by the intersection with a horizontal or vertical line respectively. 
When constraining a — 1, which gives NFW fits, the best-fit c lies in the range 2 ^ c ^ 5. 
For c = 5, the best fit a lies in the rang^j] a ^ 1. The results from the Gadget2 runs 
(GR03, G256 and G512) were in excellent agreement with the R03 runs. 

By using an analysis of the peak circular velocity, R03 found that recently virialized 
halos are well described by a constant concentration c = 5. This value was then used to break 
th e profile fitting degeneracy , giving values of a which agree with the theoretical prediction 



of ISubramanian et al.l (120001 ) for a range of halo masses from the various box sizes - for more 
details see the original R03 paper. When one uses circular velocity fitting, as in R03, with 
a Poisson assumption we obtain a smaller Act = 0.1 range. If we use density fitting with 
the NFW a = 1 constraint, illustrated by a vertical line in Figu re [p the firs t 40 h alos in 



R03-1 give c = 3.7 ± 1.4. This closely matches the results given in lColm et al.l (120041 ) - who 



analyze their lh ^^Mpc box at redshift z = 3 instead of z = 10 - once corrected for redshift 



as c oc (1 + z)^ 1 as in iBullock et al.l (120011 ). The remainder of the l/i _1 Mpc results are 
shown in Tabled! all seem to be consistent and allow for fitting using either constraint. The 
increased mean \ 2 m the G512-1 simulation occurs because the Poisson errors are reduced 
(more particles per bin) while the systematics due to departures from the fitted split power- 
law profile are constant, since these arise through "physical" processes which are not greatly 
changed in magnitude by the change in resolution. 

We now consider whether there is any evidence for differences between the l/i -1 Mpc 
and boxes and their 32 and 256/i _1 Mpc counterparts. The line of minimum x 2 f° r the 10 
most massive regular halos in R03-32 and R03-256 is plotted in the middle and right panels 
of Figure [T] respectively. Taking the 40 most massive halos and their intersections with the 
constraint c = 5 gives a = 0.6 ± 0.4 (L = 32/i" 1 Mpc) and a = 1.4 ± 0.4 (L = 256/i" 1 Mpc). 
The NFW a = 1 constraint gives c = 3.6 ± 1.4 and c = 9.1 ± 2.6 respectively. Comparing 
with the L = l/i _1 Mpc results, there is a marginal difference in the L = 32/i~ 1 Mpc results, 
and a clear difference in the L = 256/i -1 Mpc results. The comparison simulations G256-32, 
G256-256 and GR03-256 show the same trend with box size (see Tabled]). 

It is important to note that neither the NFW profile plotted in Figure 3 of R03, nor 



1 Note that in lRicotti fc Wilkinsonl (|2004r ) it is found that a concentration parameter of c = 5 is a better 
estimate than the original value c = 7 in R03. We therefore adopt this value for obtaining numerical results 
(Table [J). 
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the profiles with a shallower central slope a ~ 0.4 have reached their asymptotic values at 
the convergence radii. The innermost resolved radii, however, satisfy (3 = d In v C i r /d In r = 
1 — a/2 > 0.5, which cannot be obtained by a true NFW profile (which must have (3 < 
0.5). Unfortunately, using the density profile fitting we find that direct estimation of an 
individual profile's logarithmic slope yields results too noisy to perform direct analysis near 
the convergence radius, even with the higher resolution G512 simulations. 



3.1. A Non-parametric Method 



To circumvent the problems of the fitting degeneracy, we follow the method of lNavarro et al. 



(120041 ) . This procedure is similar to the one used by R03 where the circular velocities where 
analyzed instead of the density profiles. We make the assumption that the total mass in- 
terior to the convergence radius, r c , is representative of the "physical" halo, even if the 
inner profile is wrong. This being the case, 7 = 3 (1 — p{r)/p{r)) gives a robust limit to 
the asymptotic power law of the halos: we cannot have a steeper than p oc r -7 cusp. The 
7 - radius relation is plotted for some sample halos in Figure [2j The G512-1 halos are 
resolved to smaller radii than the R03 run. The sets are in approximate agreement where 
their data overlaps, but the innermost point of the G512 profiles have 7 < 1, suggesting 
that the profile interior to r c would not be well describ ed by NFW but probably in agree- 



ment wiht the latest high- resolution simulation results (INavarro et al.l 12004 ; iGraham et al. 



20051 ). H owever, r c is p r obabl y smaller than the core s whi ch may have b e en ob served in LSB 



see e.g. 



Kle vna et a 



ta3|), Ide Blok et all (I200lh and iGoerdt et al.l (120061 )) . That said 



Ricotti fc Wilkinson! (120041 ) have shown that the density profiles do provide good fits to the 



kinematic data from the Draco Local Group dSph. 



4. Conclusions & Discussion 



In R03, Ricotti have claimed that central slopes of dark matter density profiles (at a 
fixed fraction of the virial radius, dictated by the resolution of the simulation) are correlated 
with halo mass and that the halo concentration at the redshift of virialization is a universal 
constan t. This result is con s istent with the findin g of recent high-resolution simulations at 
z = (jNavarro et al. 2004 ; Graham et al. 2005 ) in which is found that the slope of the 
density profile in halos become flatter with decreasing radius without converging to any 
asymphotic value. 



It's often a misconception that the results in R03 contradict several published papers. 
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For example, Merritt et al. (2005, 2006) found that the average inner slope for four dwarfs 
at z = with mass 10 1 M is —1.24 as compared with —1.17 for 8 clusters. This does not 
contradict our result because clusters and dwarfs virialize at very different redshifts. The 
virial radius and concentration of dwarfs at z = has increased by a factor Zf + 1 ~ 10 from 
the redshift of their formation, Zf. Hence, the inner slope we measure, i.e., at 5%-10% of 
the virial radius in 10 8 M Q halos at z — 10, should be compared to the slope at 0.5%-0.1% 
of the virial radius at z = in comparable mass dwarfs. In addition, in the abstract of the 
same paper, they confirm a dependence of the halo profile on the halo mass, which is the 
main result in R03. Diemand, Moore, Steidel (2005), found that in minihalos of mass 10~ 6 
M at z = 26 the density profile has inner slope —1.2 at 10% of the virial radius. But the 
number of particles in the three halos they look at is a factor of 10 — 100 smaller than in our 
simulations. In addition the power spectrum of density fluctuations they use differs from 
ours as they use a power law with an exponential cut off at 0.6 parsec. It is not clear how 
the cut off in the power spectrum can affect the density profile. For example, contrary to 
our simulations, their halos do not have substructure. 

Here, we have re-analyzed the results presented in R03 and run new simulations showing 
that R03 results are reproducible by using different N-body codes, higher-resolution simula- 
tions, and different analysis techniques. Part of the tension between R03 and previous results 
may be attributed to the method of analysis of the profiles. When fitting the density profiles, 
as is done in most of previous works, the fitting degeneracy between a and c does not allow 
one to understand whether is th e slope a or the concentrat ion c that depends on the mass of 
the dark halo. If one favors the ISubramanian et al.l (120001 ) scaling arguments which express 
a relation between a and the initial power spectrum (thus, between a(r) and the enclosed 
mass M(< r)), it is found that the concentration of halos at the redshift of virialization is a 
universal constant. However, Figure [1] shows that the R03 results can be re-parameterized 
by taking a different cut through the degeneracy, leading to the more wi dely accepted notion 



of the NFW con centration depending on the halo mass and epoch (e.g lBullock et al.l 12001 



Colin et al.ll2004l ). R03, by analyzing the circular velocities of small mass galaxies at z ~ 10, 
finds some evidence for cusps shallower than a — 1 at radii < 10% of the virial radius (in 
contradiction with NFW profiles). Here, using a different non-parameterized analysis we 
also find some evidence for a flatter than a = 1 power-law of the density profile within 
5 — 10% of th e virial radius in galaxies w ith mass ~ 10 9 M at z ~ 10. Note that the scaling 
argumennt in ISubramanian et al.l (120001 ) advocated by R03 suggests that the density profile 
does not converges toward any given asymptotic value of the logarithmic slope but rather 
becomes gradually flatter toward the cente r. This is in good agreement with the re sults of 
recent high-resolution simulations at z = (INavarro et al.ll2004l ; iGraham et al.ll2005l ). How- 
ever, the present work does not have sufficient resolution to investigate this hypothesis and 
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compare the goodness of different shapes for the fitting function. 
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Table 1. 



Sim I Constraint: c = 5 Constraint: a = 1 



Mean Variance (x 2 ) Mean Variance (x 2 ) 
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3.7 
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GR03-1 
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= 0.4 
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15 
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3.8 


1.5 
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G256-1 


a 


= 0.4 
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3.4 


1.0 
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G512-1 
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= 0.2 
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32Mpc at z 
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G256-32 


a 


= 0.6 


0.4 


15 


c = 


3.8 
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256Mpc at z 


= 
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GR03-256 
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14 


G256-256 


a 


= 1.6 


0.4 


16 


c = 
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Note. — Results for fitting a split power law profile to the 40 most massive halos in each simulation, with constraints as 
follows: (1) Fix c = 5 (as in Ricotti & Wilkinson (2004)), fi = 3 and measure a; (2) Fix a = 1 (as NFW/C04), = 3 and 
measure c. There is no detectable correlation between the halo masses and any of the parameters analyzed here. Therefore 
the results are a true reflection of the ensemble of halos and their intrinsic scatter. Note that the (x 2 ) for both constraint are 
comparable: there is no evidence that one is to be preferred over the other. The higher (x 2 ) in the G512 halos is to be expected 
due to the constant systematics but smaller Poisson errors. The values obtained for the NFW concentration are in agreement 
with those obtained in other published results, once corrected for redshift (see e.g. Bullock et al. 2001) 
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Fig. 1. — From left to right: line of minimum x 2 for fits (by density) to the first 10 R03 halos 
for box sizes L — 1, 32, 256/i~ 1 Mpc. The results from the Gadget2 runs (GR03, G256 and 
G512) are in excellent agreement; G256 results are overplotted as an example of these. The 
intersection of each line with the vertical or horizontal dotted lines gives the fit constrained 
to an NFW profile or a fixed c split power-law profile respectively. While x 2 varies along each 
line, it does so by a factor of only a few; moreover the mean x 2 obtained for constraining a 
number of halo fits along one axis or the other is almost identical (see Table [1]). 




100 1000 
r/pc 

Fig. 2. — The first five most massive G512 (dashed) and R03 halo 7 = 3(1 — p/p) profiles. 
Each profile is plotted only exterior to its convergence radius (which are smaller for the G512 
simulations). 7 would converge to a in the case of the split power law profile. There does 
not appear to be any evidence for convergence to a particular value. 



